{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5678b63a-3352-4256-8691-d9578a292fa9",
   "metadata": {},
   "source": [
    "## Decomposition after structural estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "56a0e031-ed12-4bb1-9820-7183d7a4c0e4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from fn_decompose import *\n",
    "\n",
    "# randomization variables for simulation\n",
    "#import secrets\n",
    "#secrets.randbits(128) \n",
    "seed = 84708628852577542122415832887111152745\n",
    "\n",
    "rng = np.random.default_rng(seed)\n",
    "\n",
    "n_draws = 1000\n",
    "\n",
    "param_random = {'n_draws' : n_draws, \n",
    "                'seed' : seed}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "91552552-8828-48c0-88a2-3c69a7fd5875",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# calibrated parameters\n",
    "\n",
    "    # income\n",
    "mu_lny = 9.713\n",
    "sig_lny = 0.497\n",
    "\n",
    "mu_eps = 0\n",
    "\n",
    "    # prices\n",
    "pi_s = 1\n",
    "pi_q = 0\n",
    "tau = 0.075\n",
    "\n",
    "gamma = 0.424\n",
    "\n",
    "param_given = {'mu_lny' : mu_lny, \n",
    "               'sig_lny' : sig_lny,\n",
    "               'mu_eps': mu_eps,\n",
    "               'pi_s' : pi_s, \n",
    "               'pi_q' : pi_q, \n",
    "               'tau' : tau,\n",
    "               'gamma': gamma}\n",
    "\n",
    "# estimated parameters\n",
    "# initial: [ 1.467219  0.110516  0.92719   0.611335 -7.924372  0.299971 -0.139612 -0.010366]\n",
    "# results:  [ 1.466967  0.110715  0.927835  0.612176 -7.919992  0.299912 -0.140181 -0.009912]\n",
    "\n",
    "# results:\n",
    "pi_n = 1.466967\n",
    "pi_nq = 0.110715\n",
    "\n",
    "theta = 0.927835\n",
    "alpha = 0.612176\n",
    "rho = -7.919992\n",
    "\n",
    "sig_eps = 0.299912\n",
    "eps_cutoff = -0.140181\n",
    "\n",
    "cor_lny_eps = -0.009912\n",
    "\n",
    "\n",
    "param_choice = np.array( \n",
    "    [pi_n, pi_nq, \n",
    "     theta, alpha, rho,\n",
    "     sig_eps, eps_cutoff, cor_lny_eps] \n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "71aa1bb0-de2f-4a9c-837e-09d71d214154",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "disp = False\n",
    "\n",
    "# randomization variables\n",
    "n_draws = param_random['n_draws']\n",
    "seed = param_random['seed'] # seed = a integer or None\n",
    "\n",
    "mu = [mu_lny, mu_eps]\n",
    "cov_lny_eps = cor_lny_eps * sig_lny * sig_eps\n",
    "\n",
    "cov = [[sig_lny**2,cov_lny_eps],[cov_lny_eps,sig_eps**2]]\n",
    "\n",
    "    # draw from a joint-normal distribution of log(income) and eps\n",
    "rng = np.random.default_rng(seed) # use a pre-specified seed for replication\n",
    "lny_eps_draws = rng.multivariate_normal(mu, cov, n_draws)\n",
    "\n",
    "# container of decomposition results\n",
    "dData_list = []\n",
    "\n",
    "# create an instance of FamilyRation()\n",
    "FR = FamilyRation(\n",
    "    pi_n = pi_n, pi_nq = pi_nq, pi_s = pi_s, pi_q = pi_q, \n",
    "    theta = theta, rho = rho, gamma = gamma, \n",
    "    alpha = alpha, y = 10\n",
    "    )    \n",
    "\n",
    "# main simulation, repeated by \"n_draws\" times\n",
    "if disp: \n",
    "    print(\"Running a simulation ... \" + str(n_draws) + \" iterations\")\n",
    "for i in range(n_draws):\n",
    "\n",
    "    [lny, eps] = lny_eps_draws[i]\n",
    "\n",
    "    if disp: \n",
    "        print(\"---------------------\")\n",
    "        print(\"Iteration: \", i)\n",
    "        print(\"[lny, eps]: \", lny, eps)\n",
    "        print(\"---------------------\")\n",
    "\n",
    "    # income exponential and scaling\n",
    "    y = np.exp(lny)/1000\n",
    "\n",
    "    # add a proportion of income to child price\n",
    "    pi_n_tau = pi_n + tau * y\n",
    "\n",
    "    # check income\n",
    "    if y - pi_n_tau * 3 < 0.5:\n",
    "        if disp: print(\"Insufficient income to bear three children, pass\")\n",
    "        continue\n",
    "\n",
    "    # adjust preference parameter alpha\n",
    "    theta_new = theta + eps\n",
    "    if theta_new < 0:\n",
    "        theta_new = 0.01\n",
    "    elif theta_new > 1:\n",
    "        theta_new = 0.99\n",
    "\n",
    "    # update FR\n",
    "\n",
    "        # symbolic\n",
    "    FR.y = y\n",
    "    FR.pi = np.array( [pi_n_tau, pi_q, pi_s, pi_nq] )\n",
    "    FR.a = np.array([alpha, rho, theta_new, gamma])\n",
    "\n",
    "    # decompose\n",
    "    dData = decompose(FR, eps, eps_cutoff, disp = disp)\n",
    "    \n",
    "    dData_list.append(dData)\n",
    "\n",
    "# generate simulation dataframe\n",
    "dData_df = pd.concat(dData_list).reset_index(drop=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c808df72-ea7b-45f6-bd9f-d5a8d427c18f",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>type2to3</th>\n",
       "      <th>q_2half</th>\n",
       "      <th>total_2half</th>\n",
       "      <th>price_2half</th>\n",
       "      <th>substitution_2half</th>\n",
       "      <th>income_2half</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>C</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>C</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>C</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>A</td>\n",
       "      <td>8.196025</td>\n",
       "      <td>0.242021</td>\n",
       "      <td>-0.121562</td>\n",
       "      <td>0.432067</td>\n",
       "      <td>-0.068580</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>C</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>967</th>\n",
       "      <td>C</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>968</th>\n",
       "      <td>B</td>\n",
       "      <td>10.696154</td>\n",
       "      <td>0.446661</td>\n",
       "      <td>-0.105932</td>\n",
       "      <td>0.355775</td>\n",
       "      <td>0.196214</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>969</th>\n",
       "      <td>A</td>\n",
       "      <td>7.509682</td>\n",
       "      <td>0.277083</td>\n",
       "      <td>-0.139841</td>\n",
       "      <td>0.571573</td>\n",
       "      <td>-0.154968</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>970</th>\n",
       "      <td>A</td>\n",
       "      <td>2.882801</td>\n",
       "      <td>-1.242381</td>\n",
       "      <td>-0.221391</td>\n",
       "      <td>0.124705</td>\n",
       "      <td>-1.145660</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>971</th>\n",
       "      <td>C</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>972 rows × 6 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "    type2to3    q_2half  total_2half  price_2half  substitution_2half  \\\n",
       "0          C        NaN          NaN          NaN                 NaN   \n",
       "1          C        NaN          NaN          NaN                 NaN   \n",
       "2          C        NaN          NaN          NaN                 NaN   \n",
       "3          A   8.196025     0.242021    -0.121562            0.432067   \n",
       "4          C        NaN          NaN          NaN                 NaN   \n",
       "..       ...        ...          ...          ...                 ...   \n",
       "967        C        NaN          NaN          NaN                 NaN   \n",
       "968        B  10.696154     0.446661    -0.105932            0.355775   \n",
       "969        A   7.509682     0.277083    -0.139841            0.571573   \n",
       "970        A   2.882801    -1.242381    -0.221391            0.124705   \n",
       "971        C        NaN          NaN          NaN                 NaN   \n",
       "\n",
       "     income_2half  \n",
       "0             NaN  \n",
       "1             NaN  \n",
       "2             NaN  \n",
       "3       -0.068580  \n",
       "4             NaN  \n",
       "..            ...  \n",
       "967           NaN  \n",
       "968      0.196214  \n",
       "969     -0.154968  \n",
       "970     -1.145660  \n",
       "971           NaN  \n",
       "\n",
       "[972 rows x 6 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dData_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2f94e9c8-216c-429d-8430-d6d487906c9e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "pd.concat(dData_list).to_csv('./decompose_main.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa294480-4a36-4af4-b888-bfb0a7bb409a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:base] *",
   "language": "python",
   "name": "conda-base-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
